!***********************************************************************
! Equation Solving Module
!
! 20/Feb/2013
! Functions -
!   1. phyGaussJordan(matrix, solution)
! 2. phyGaussElemination(matrix, solution)
!   3. phyCheckError(matrix, solution, mean)
!   
!-------------------------------------------------------------------    
! 
! 1. Gauss - Jordan Equation Solver
!
! Finds the solution to the matrix  by diagonalizing the matrix of coefficients
!   Uses pivoting to reduce the error due to truncation of floating points
!   usage - call phyGaussJordan(matrix, solution)    
!   matrix   -> (n,n+1) matrix where n-rows, n+1 columns, the nxn matrix should contain the coefficients
!              the n+1 column should contain the constants   
! solution -> A 1xn column vector, which will contain the solution of the system of equations. Incase no solution
!               is found. It will be set to zero. Call the function phyErrorCode() to see what error has occured
!----------------------------------------------------------------------------------------------------------------
!
! 2. Gauss - Elimination Equation Solver
!
! Finds the solution to the matrix  by upper triangularizing the matrix of coefficients
!   Uses pivoting to reduce the error due to truncation of floating points
!   Usage - call phyGaussElemination(matrix, solution)   
!   matrix   -> (n,n+1) matrix where n-rows, n+1 columns, the nxn matrix should contain the coefficients
!              the n+1 column should contain the constants   
! solution -> A 1xn column vector, which will contain the solution of the system of equations. Incase no solution
!               is found. It will be set to zero. Call the function phyErrorCode() to see what error has occured
!----------------------------------------------------------------------------------------------------------------
! 3. Jacobi-Iterations
!
! Finds solutionsof a set of equations using the Jacobi Iteraton Method. The iteration stops when the mean error
! calculated by phyCheckError() reaches limit
!
!   Usage - call phyJacobiIteration(matrix, solution,limit)
!   matrix   -> (n,n+1) matrix where n-rows, n+1 columns, the nxn matrix should contain the coefficients
!              the n+1 column should contain the constants   
! solution -> A 1xn column vector, which will contain the solution of the system of equations.
! limit -> the bounds for error
!----------------------------------------------------------------------------------------------------------------
! 4. Mean Error due to Floating point truncations
!
! Finds the average error in the solution of the equations
!
!   Usage - call phyCheckError(matrix, solution,mean)
!   matrix   -> (n,n+1) matrix where n-rows, n+1 columns, the nxn matrix should contain the coefficients
!              the n+1 column should contain the constants   
! solution -> A 1xn column vector, which will contain the solution of the system of equations.
! mean       -> Returns the mean error.typically of the order of 10^-3
!----------------------------------------------------------------------------------------------------------------
! SAMPLE CODE
! program test
!     use EquationSolving_module
!     implicit none    
!     integer             ::  i, j, N
!     real, allocatable   ::  matrix(:,:), sol(:)
!     real                ::  meanerr

!     print*,"Enter the number of variables"; read*,N
!     allocate(matrix(N,N+1))
!     allocate(sol(N))
!     print *, 'enter the matrix';    
!     do i = 1, N 
!         read *, (matrix(i,j),j=1,N+1)
!     end do
!     call phyGaussElemination(matrix, sol)
!     print*,sol
!     call phyCheckError(matrix, sol, meanerr)
!     print*,"Mean Error=",meanerr
! end program
!===================================================================================================================

module EquationSolving_module
    implicit none
    integer :: EquationSolving_error
contains
    
    subroutine  phyGaussJordan(matrix, solution)
        real, allocatable, dimension(:,:), intent(in)   ::  matrix
        real, allocatable, dimension(:), intent(out)    ::  solution
        real, allocatable, dimension(:,:)               ::  gauss_mat
        integer     ::  rows, cols
        integer     ::  i, j, index_of_pivot
        integer     ::  m,n
        EquationSolving_error = 0
        rows = size(matrix, dim=1)
        cols = size(matrix, dim=2)
        allocate(gauss_mat(rows, cols))
        gauss_mat = matrix
        !generating the diagonal form selecting the pivot
        do j = 1, cols-1
            index_of_pivot = absMaxOf(gauss_mat(j:rows,j))  !pivot is the absolute maximum of the jth column            
            index_of_pivot = index_of_pivot + j - 1
            call swapRows(gauss_mat, j, index_of_pivot)
            if( gauss_mat(j,j) == 0 ) then
                allocate(solution(rows))
                solution = (/(0, i=1,rows)/)
                !print*,"No solution"
                EquationSolving_error = 1               
                return
            end if      
            gauss_mat(j,:) = gauss_mat(j,:)/gauss_mat(j,j)
            do i = 1, rows
                if(i == j) cycle
                gauss_mat(i,:) = gauss_mat(i,:) - gauss_mat(i,j)*gauss_mat(j,:)
            end do
        end do
        allocate(solution(rows))
        solution = gauss_mat(:,cols)
        EquationSolving_error = 0
    end subroutine phyGaussJordan

subroutine  phyGaussElemination(matrix, solution)
        real, allocatable, dimension(:,:), intent(in)   ::  matrix
        real, allocatable, dimension(:), intent(out)    ::  solution
        real, allocatable, dimension(:,:)               :: gauss_mat
        integer     ::  rows, cols
        integer     ::  i, j, index_of_pivot, m
        rows = size(matrix, dim=1)
        cols = size(matrix, dim=2)
        allocate(gauss_mat(rows, cols))
        gauss_mat = matrix
        
        !generating the upper triangular form selecting the pivot   
        do j = 1, cols-1
            index_of_pivot = absMaxOf(gauss_mat(j:rows,j))  !pivot is the absolute maximum of the jth column            
            index_of_pivot = index_of_pivot + j - 1
        do m = 1, rows
      end do
            call swapRows(gauss_mat, j, index_of_pivot)
            if( gauss_mat(j,j) == 0 ) then
                allocate(solution(rows))
                solution = (/(0, i=1,rows)/)
                !print*,"No solution"
                EquationSolving_error = 1
                return  
            end if      
            gauss_mat(j,:) = gauss_mat(j,:)/gauss_mat(j,j)
            do i = j+1, rows
                gauss_mat(i,:) = gauss_mat(i,:) - gauss_mat(i,j)*gauss_mat(j,:)
            end do
        end do
        allocate(solution(rows))
        solution = (/(0,i = 1,rows)/)
        solution(rows) = gauss_mat(rows,cols)
        do i=rows-1, 1, -1
            solution(i) = gauss_mat(i,cols)-DOT_PRODUCT(gauss_mat(i,:),solution(:))
        end do
        EquationSolving_error = 0
    end subroutine phyGaussElemination

    subroutine phyJacobiIteration(matrix, solution, limit)
        real, allocatable, dimension(:,:), intent(in)   :: matrix
        real, allocatable, dimension(:), intent(out)    :: solution
        real, allocatable, dimension(:,:)               :: b, temp, temp2        
        real, allocatable, dimension(:,:)               :: D, R
				real, intent(in)																:: limit        
				integer :: i, j
        integer :: N
        real    :: err
        N = size(matrix, dim=1)
        allocate(D(N,N))
        allocate(R(N,N))
        allocate(b(N,1))
        allocate(temp(N,1))
        allocate(temp2(N,1))
        allocate(solution(N))
        do i = 1,N
            solution(i) = 0
            temp(i,1) = 0
            temp2(i,1) = 0
            D(i,:) = (/(0,j=1,N)/)
            D(i,i) = matrix(i,i)
            R(i,:) = (/(0,j=1,N)/)
        end do
        R = matrix(:,:N) - D
        b(:,1) = matrix(:,N+1)
        !taking the inverse of the diagonal matrix
        do i = 1, N
            D(i,i) = 1/matrix(i,i)
        end do
        !applying the iteration X(k+1) = D(^-1)(b - RX(k))
        do    
            temp2 = MATMUL(D,b  - MATMUL(R, temp))
            solution = temp2(:,1)
            call phyCheckError(matrix, solution, err)
            if (abs(err) < limit) exit
            temp = temp2
        end do
        solution = temp2(:,1)
    end subroutine phyJacobiIteration

    subroutine phyCheckError(matrix, solution, mean)
        real, allocatable, dimension(:,:), intent(in)   :: matrix
        real, allocatable, dimension(:), intent(in)     :: solution
        real, intent(out)                               :: mean             
        real, allocatable, dimension(:)                 :: error
        integer :: i, j
        integer :: rows, cols
        real        :: sum
        rows = size(matrix, dim = 1)
        cols = size(matrix, dim = 2)
        allocate(error(rows))
        sum = 0
        do i = 1,rows
            error(i) = abs(DOT_PRODUCT(matrix(i,:cols-1),solution) - matrix(i,cols))
            sum = sum + error(i)
        end do
        mean = sum/real(rows)
    end subroutine phyCheckError

    subroutine  swapRows(matrix, m, n)
        real, allocatable, dimension(:,:), intent(inout)    ::  matrix
        integer, intent(in)                             ::  m, n
        real, allocatable, dimension(:)     ::  temp
        allocate(temp(size(matrix, dim=2)))
        temp        = matrix(m,:)
        matrix(m,:) = matrix(n,:)
        matrix(n,:) = temp
    end subroutine swapRows

    function absMaxOf(column)
        real, dimension(:), intent(in)  ::  column
        real        ::  temp(size(column))
        integer ::  absMaxOf
        integer ::  i
        real        ::  max
        temp = abs(column)
        absMaxOf = maxloc(temp, 1)
    end function absMaxOf
end module  EquationSolving_module
